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We study the decay of multiple quantum (MQ) NMR coherences in systems with the large number 
of equivalent spins. As being created on the preparation period of MQ NMR experiment, they decay 
due to the dipole-dipole interactions (DDI) on the evolution period of this experiment. It is shown 
that the relaxation time decreases with the increase in MQ coherence order (according to the known 
' results) and in the number of spins. We also consider the modified preparation period of MQ NMR 

^ ■ experiment (G.A.Alvarez, D.Suter, PRL 104, 230403 (2010)) concatenating the short evolution 

■ periods under the secular DDI Hamiltonian (the perturbation) with the evolution period under the 



non-secular averaged two-spin/two-quantum Hamiltonian. The influence of the perturbation on the 
decoherence rate is investigated for the systems consisting of 200-600 equivalent spins. 
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I. INTRODUCTION 

^ ^, Multiple quantum (MQ) coherences are quite suitable for investigations of the dependence of the relaxation time 



on the size of the quantum system This problem is closely connected to the estimations of the decoherence 

' time as an important parameter for the quantum information systems. A simplest model of the quantum register 

> ■ . n 

' formed by the highly correlated spins can be created in MQ NMR experiments [5[- Some models of quantum 

OS 



The theoretical methods (for example, 



registers consisting of up to 4900 qubits were studied experimentally 
I ref.Q) describing the experiments 1| are phenomenological ones and the development of theoretical and numerical 

00 

. allow us, generally speaking, to investigate systems consisting of not more than twenty spins [7|. Some progress in 
the study of the larger systems (up to 40 spins) is achieved due to the special techniques based on the Chebyshev 
• polynomial expansion 0, Q] and on the phenomenon of quantum parallelism [lol | . The new pers pec tives arc opened 
' by MQ NMR in systems of equivalent spins where the special method has been worked out [ill. I12I allowing one to 
investigate MQ NMR dynamics of hundreds of spins and even more. Such systems of equivalent spins can be created 
in a nanopore compotirid placed in a strong external magnetic field if the nanopores are filled with a gas of spin-carring 
IJ] . Since the characteristic time of the molecular diffusion is much less than the spin flip-flop 
14l |. the dipole-dipole interactions (DDI) of spins arc averaged ( but not to zero) and the residual DDI can 
be described by the single coupling constant. As a result, all spins can be considered as equivalent ones, which 
significantly simplifies the numerical simulations. 



molecules (atoms) 



time 



13 



'Electronic address: 's . i . doronin@gm ail.com| 
tElectronic address: cfeldman@icp. ac.ru| 
t Electronic address: |zenchuE@itp.ac.ru| 



2 



H 



eff 



H 



ev 



preparation, l evolution, t 



mixing, r 



detection 



FIG. 1: The basic scheme of the MQ NMR experiment. The Hamihonian Hefs (Eq.©), Hev (Eq.(|8l))and Hmq (Eq.©) govern 
the spin dynamics on the appropriate period of the MQ NMR experiment 



The above method can be applied to the investigation of the decay of MQ NMR coherence intensities of different 
orders caused by the secular DDI in systems containing hundreds of spins. In the simplest case this decay occurs 
on the evolution period of the MQ NMR experiment However, the MQ NMR experiment can be modified, for 
instance, using a different set of pulses on the preparation period, as in Ref.Q. In the later case, the decay takes 
place on the preparation period. 

In this paper we study the decay of MQ NMR coherence intensities created in systems with a large number of 
equivalent spins. The paper is organized as follows. The general description of different MQ NMR experiments is 
given in SecHIl The theory and the numerical simulation of the decay of MQ NMR coherence intensities in different 
MQ NMR experiments is developed in Scc lIIII The conservation law associated with considered models is derived in 
Sec lIVI We briefly summarize our results in concluding SecfVl 



II. THE MQ NMR EXPERIMENTS IN A SYSTEM OF EQUIVALENT SPINS 

The MQ NMR experiment consists of four distinct periods of time (FiglT]): preparation (r), evolution {t), mixing 
(r) and detection. 

a. Preparation period. The spin system is irradiated by the proper multipulse sequence on the preparation period. 
As a result, the anisotropic DDI of nuclear spins in the external magnetic field, i?, (directed along the axis z) oscillates 
rapidly. In the rotating reference frame [15| , the dynamics of spin system is described by the effective Hamiltonian 
Hef f ■ We consider two types of pulse sequences on the preparation period. The first one is the standard pulse sequence 
resulting in the averaged non-secular two-spin/two-quantum Hamiltonian, describing the MQ NMR dynamics on the 
preparation period of the standard MQ NMR experiment , i.e. Heff = Hmq'- 

Hmq=H^+'^+H^-^\ (1) 

Here Djk = 7^^(1 ~ 3 cos^ ^^jfe)/(2r|J,) is the coupling constant between spins j and k, 7 is the gyromagnetic ratio, 
rjk is the distance between spins j and k, 6jk is the angle between the vectors fjk and B and ~ Ij^ ± iljy are 
the raising and lowering operators of spin j. The second type of pulse sequences is introduced in Ref.Q, where a 
modification of the preparation period of the MQ NMR experiment was suggested. In this case, the preparation 
period consists of the cycles of the duration Tc and each cycle concatenates the short evolution period Tdz under the 
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perturbation Hamiltonian Hdz (which is responsible for the secular DDI 15|), 



Hdz — ^ DjkC^Ijzhz — Ijxhx — Ijyhy) (3) 

j<k 

with the evolution period tmq under the ideal MQ Hamiltonian Hmq (IIJ- Thus, Tc = Tdz + tmq. Introducing the 
relative strength p = Tdz/Tc (0 < p < 1) of the perturbation one can find that the resulting evolution can be described 
by the effective Hamiltonian H^f / given by the following equation : 

H,ff {p) = {l- p)Hmq + pHdz. (4) 

Let the preparation period consist of K (a big number) cycles of the duration Tc, so that one can introduce the 
parameter r = Kt,,. Note, that Hgff{0) = Hf^jQ, which means that the standard preparation period, used, for 
instance, in ref. [^, is a particular case of the described modification. 

Hereafter we study the MQ NMR dynamics of equivalent spins. Such a case can be realized, for instance, by the 
dipolar c oupling spins in a nanopore where the Hamiltonian ([l} is averaged ( but not to zero) by the fast molecular 



diffusion 
as follows 
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14 L The Hamiltonians Hmq and Hdz with the averaged coupling constant D {Dtc <C 1) can be rewritten 
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Hmq = -j{{n' + (5) 
H,z^^{3l'z-I'}, (6) 

N N 

where ~ if [N is the number of spins), Iz ^ Y ^zj ^^'^ operator P is the square of the total spin angular 
momentum. Thus Eq.(|3]) must be replaced with the following one 

H,ff{p) = {1-p)Hmq +pHdz, (7) 

which is valid for the system of equivalent spins. 

b. Evolution period. Let the spin system be governed by the following general Hamiltonian H^y on the evolution 
period: 

H,, = {l- 9{p))Hdz + AIz, e={ . (8) 

1, p>0 



The offset A encodes MQ NMR coherences of different orders, see below, Eq. p4|) and ref.[ 

c. Mixing period. The spin system on the mixing period is governed by the Hamiltonian —Hmq in all experiments 
considered in this paper. 

We emphasize that the decay of MQ NMR coherences is caused by the Hamiltonian Hdz appearing either on the 
evolution period (if p = in Eqs. ()7|8p ) or on the preparation period (if p 7^ in eqs. ()7|8|) ). These two cases are 
considered separately in Secs lIII A\ and HITB] respectively. 
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III. THE DECAY OF MQ NMR COHERENCE INTENSITIES CAUSED BY THE SECULAR DDI 
A. The decay of MQ NMR coherences in MQ NMR experiments of Ref. Ql 

We consider the time evolution of the coherences in MQ NMR experiments with the standard preparation period 
For this purpose we take p — in Eqs.^ and dH), which read: 

Heff = Hmq, Hev = Hdz + AI,, (9) 
so that the coherence decay occurs on the evolution period. In order to investigate the MQ NMR dynamics of the 



system one should find the density matrix p{t) on the preparation period solving the Liouville evolution equation 15 1 

^^^iHMQ,p(.r)] (10) 



with the initial thermodynamic equilibrium state p{0) = Iz in the high temperature approximation 15|. Taking into 
account the pointed information about the Hamiltonians on the different periods of MQ NMR experiment one can 
write the expression for the longitudinal polarization (Iz) (t, t) after the mixing period of MQ NMR experiment (FigH]) 
as follows: 

{Iz){T,t) = Tr{C/+(T)e-^^*^^e-'^-*C/(r) x (11) 
4C/+(T)e*^"-e'-f^'^-*C/(T)/J = 

where p(t) = U{t)IzU^{t) is the solution to Eq. ([TO| and U{t) = exp{~iHMQT). It is convenient to expand the spin 
density matrix p{t) in the series as follows 

p{r)=Y.P^(^)' (12) 

fc 

where PfejV) is the contribution to p{t) from MQ coherence of the fcth order and satisfies the following commutation 



relation 



III: 



e pkc — e pk- [i-6) 

Then Eq.([ni) reads 

{Iz){T,t) = ^e-''=^*Tr{e-'^^-Vfc(r)e^^-V-fc(T)}. (14) 
fc 

Eq. ([T4| defines the intensity Jfc(T, t) of the MQ NMR coherence of order k as follows: 

Jfc(r, t) = Tr{e-'^- Vfe We'^^^'^V-fcM}. (15) 



In analogy to the autocorrelation function for the decay of the transverse magnetization [15[, Eq.(|T5|) reveals the decay 
of MQ NMR coherences due to the secular DDI on the evolution period. Since we consider a system of equivalent 
spins, numerical simulation of Eq. (jl5p may be simplified allowing one to perform calculations in the systems with the 
large number of spins. This happens due to the commutation relation 

[Hmq,I^]^0, (16) 



llj . It was shown 



that the Hamiltonians 



which suggests us to use the basis of common eigenvectors of and 
Hmq, Hdz and the density matrix p(t) for the system of equivalent spins have a block structure. For instance, 
Hmq = dia.g{H^jQ, H^^q , ■ ■ ■ , H^jq ^ }, ([a] is an integer part of a) . These blocks correspond to different total 
spin numbers S = N/2, N 12 — 1, . . . , — [N/2] [l^. All blocks are degenerated and their degeneracy nN{S) is 



determined as follows 
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^•!(25+l) ^ ^ N 



(17) 



2 ' ^ ' 2 

Thus, the problem is reduced to the set of analogous problems of lower dimensions. The intensities of MQ NMR 
coherences Jk,s{T,t) can be calculated for all blocks. Then the observable intensities Jk{T,t) {—N < k < N) are 
following 



Hi: 



Jk{T,t) = ^nN{S)Jk,s{T,t). 



(18) 



The results of numerical simulations of these intensities are represented below. 



1. The numerical simulations 

We study the dynamics of MQ NMR coherence intensities in the nanopore filled with the spin-carring particles 
numerically. Let us emphasize one more time that wc arc dealing with the highly symmetrical model where any two 
spins interact with the same constant of DDI because the diffusion characteristic time in the nanopore is much shorter 
then the spin flip-flop time [l^ Ql ■ This fact simplifies the numerical calculations significantly since all particles are 
"nearest neighbors" in this model and we consider interactions among all of them. 

Our calculations showed [lll l th at MQ NMR coherence intensities are quickly oscillating functions. For this reason 

and consider the averaged intensities 



we follow the strategy of Ref. 
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^ro+2T 

Mt)^l^ Jfe(r,t)df, (19) 



To 



= 31, T = 27r/A,„i„ = 27r/V3. 



where r = Dt and t ~ Dt arc the dimensionless times associated with the preparation and evolution periods 

_ 3 

respectively, Xmm is the minimal eigenvalue of the Hmq Hamiltonian. This value belongs to the block H^jq of 



the Hamiltonian |ll| . The choice of tq is motivated by the requirement that the coherences of all possible orders 
have appeared and one can think that the quasi-stationary distribution of the intensities is realized, which has been 
verified in ref.[ll|. The averaging is performed over two maximal periods T of oscilations, which is taken from the 
requirement that the increase of the averaging interval does not change The averaged intensities decay 

with the time t of the evolution period. The time moments (such that J/c(0)/Jfe(t)|t=t^ = e for an arbitrary k) 
versus MQ coherence order in systems with 201, 401 and 601 spins are shown in FigO We can see from this figure 
that 



1. MQ NMR coherence decay times decrease with the increase in the number of spins; 

2. MQ NMR coherence decay times decrease with the increase in their order. 
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FIG. 2: The decay time as a function of the coherence order for the spin systems with A'^ = 201, 401 and 601. 

The times te(fc) of the decay of MQ NMR coherences of order fc > can be approximated by the hypcrboHc cotangent, 
as it is shown in Figl^l 

te(fc) = fli coth(a2fc + as), (20) 
where parameters 01,02,03 may be found by the least square method: 



te{k) 



The approximation given by Eq. (|20|) shows that the decay rate of the high order coherence intensities is ahnost 
independent on their order. This happens because MQ coherence phases (acquired during the evohition period) are 
approximately proportional to their order, see Eq. (fT3|). As a result, the rates of MQ coherence decays increase with 
their order and the decay time te{k) is ^ 1/k. Thus, for the high order coherences, we have te{k)/te{k + 1) ^ 1, i.e. 
the decay times of the fcth and (fc + l)th coherences are almost the same, which is reflected in Eq. ((20|) . Regarding the 
zero-order coherence, its intensity Jo does not decay owing to the commutation relation [Hdz,Pa] — 0, which follows 
from the fact that both Hdz © and po are diagonal in the chosen basis. 

It is worthwile to note that the dynamics of the multi-spin cluster growth during the evolution of the solid spin 
system considered, for instance, in P, [l^ is essentially different in comparison with that in the system of equivalent 
spins. The matter is that only the strongly interacting spins are joined in the clusters initially, usually the nearest 
neighbors in the crystal lattice [sl, [l^. After that, the next neighbors become involved in the cluster and so on. Thus 
more and more remote spins become embedded in the cluster with time. As a result, it becomes possible to observe 
the growth of the multi-spin clusters in MQ NMR experiments [l^. However, the dynamics of the spin clusters 
is quite different in the high symmetrical spin system such as the system of equivalent spins. All spins are "nearest 
neighbors" in this case, so that the spin cluster consisting of all N spins is formed much more quickly during the time 
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FIG. 3: The evolution of the "cluster size" Nc\ t is the dimensionless evolution time 



interval ^ l/D, where D is the constant of DDI, which is the same for any two spins. It becomes hard to follow the 
process of the cluster growth in the high symmetrical system of equivalent spins, unlike the solids jsjjll^. Nevertheless, 
there is some reorganization of the spin cluster diring the evolution, when the system is irradiated by the multipulse 
sequence 0, Q, resulting to the high order MQ coherences. 

Now let us turn to the decay of MQ NMR coherences. We consider the "cluster" of MQ NMR coherences as a family 
of such coherences whose intensities exceed some fixed value Jmm, say, Jmm = 0.005. This minimal value is taken 
since the smaller intensities are hardly observable in the experiment. The size Nc of the cluster of MQ coherences does 
evolve, which is demonstrated in FigIS] This evolution is a consequence of the fact that the rate of the decay increases 
with the increase in the order of MQ NMR coherences. We see also that the rate of decrease of the coherence cluster 
size increases with the increase in A^. The described experiment may be used in order to prepare the coherence 
clusters of desirable size varying the duration of the evolution period. 

B. The decay of MQ NMR coherences in MQ NMR experiments with the modified preparation period 

It is very important to investigate the degradation of quantum superposition states. MQ NMR experiments jsj allow 
us to make it. To this end the modification of the preparation period of the MQ NMR experiment was suggested 
in rcf.Q. In this section we consider Eqs.Q and ([8]) with p > 0, so that the system is governed by the general 
Hamiltonian i?e// during the preparation period and by the Hamiltonian H^v = during the evolution period. 
Thus the coherence decay occurs on the preparation period of the MQ NMR experiment. The calculations analogous 
to those used for the derivation of Eq. p^ yield the following expression for the intensities of MQ NMR coherences: 



Jk{T,p) = Tr{pfe(r,p)p_fe(r,p)}, r = Ktc 



Cj 



(21) 



where 




(22) 



k 
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If p ^ 1, then it is simple to demonstrate that the intensities vary proportionally to . In fact, the Liouville 
equation on the preparation period can be rewritten as follows: 

" '^^^-p)HMQ+pHd.rp{T)\. (23) 

Solving Eq. ((23|) by the methods of the perturbation theory [l^ one can obtain 

p[t) = p{t) - ppi (r) - p^p2 (r) , (24) 
Jo 

HMQ)e-'"'"^^^'-^\p{T)]dT\ 

P2(r)= r dr' r dr" [e^^-« (ff,,- 
Jo Jo 

i^MQ)e-^^-«(^"-^\p(T)]], 

where p(r) is the solution to Eq. pi?)) . It is evident from Eg. ip^ that 

Tr{p(r)p(r)} = Tr{p'{r)} - p'A{t), (25) 
A{t) ^ J dr' J rfr"Tr| 
[p(t),e'^"«(^'-^)(i?rf, -i7MQ)e-'^"«("'-")] X 

[g,5„«(r"-r)(^^^ - i?A/Q)e-'^"«(""-^',p(r)]}. 

The behavior of the intensities with the increase in p is defined by the sign of A{t). Wc do not determine this sign 
for an arbitrary t. However, one has for small t: 

A « -yTr[p, Ha, - Hmq]^ > 0. (26) 

Since Tr{p^} is the sum of the intensities of MQ NMR coherences for the standard MQ NMR experiment and 
Tv{p{t)p{t)} is the analogous sum, when the perturbations are taken into account one can conclude that the intensities 
decrease with the increase in the square of the perturbation strength at least for small r, which is confirmed below 
by the numerical simulations. 



1. The numerical simulations 

We refer to f = Dt as the dimensionless evolution time in this section. Before proceed to the numerical simulations, 
let us underline the basic difference between the experiments considered in Secs lIII Al and HlTBl The matter is that the 
high frequency oscillations of MQ NMR coherences are formed on the preparation period with the duration r, while 
the decay of these coherences occurs on the evolution period (with the duration t) in Sec lIII Al For this reason we 
consider the intensities averaged over the parameter r in that section. However, the situation is different in Sec lIII Bl 
because the decay occurs on the preparation period, so that the parameter r is responsible for the both oscilations 
and decay of MQ NMR coherences. Because of this fact, we are not able to consider the averaged intensities. 
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Coherence order, k 

FIG. 4: The decay time as a function of the coherence number for spin systems with = 201 

Instead of this, wc relate the decay time Tp(fc) of the intensity Jk{f) with the decay times T^{k) of its envelopes. 
Here the subscript p indicates that the parameter Tp depends on the value of p in the Hamiltonian H^ff. Parameters 
T^{k) may be found simply by plotting the graphs of the envelopes J^{t) of the quickly oscillating intensity Jfe(f), 
(f) < Jfc(f) < J^(f). Then the decay times of the envelopes T^{k) are the first zeros of J^{f) appeared after 
the amplitude of the fcth intensity gets its maximal value: J^{Tp™{k)) ~ 0. All this suggests us to calculate the 
decay time Tp(fc) of the fcth coherence as follows. First, we have to find numerically all solutions Tp\k) (i = 1, 2, . . . ) 
to the equation Jfc(T,p)|__^(i)|,^^ = 0, such that Tp{k) < Tp\k) < Tj^{k). Let Np{k) be the number of such 
solutions. Then the decay time of the fcth coherence intensity may be found as the averaged value of these solutions: 
Tp{k) = jfj/Tj X^fcTi'^'' '''p\^)- The dependence of Tp(fc) on the coherence number is shown in Fig|3]for N = 201. It is 
found that this decay may be approximated as follows: Tp(fc) fa Up + bp tanh(rfp — Cpk). Parameters Op, fop, Cp, dp have 
been found by the least square method, see Table I. We see that the decay time of the high order coherence intensity 



p 


Up 


bp 


Cp 


dp 


0.001 


42.0073 


19.7734 


0.0565 


1.5240 


0.002 


21.1130 


10.5523 


0.0543 


1.4369 


0.003 


14.9127 


7.4843 


0.0472 


1.1474 


0.004 


10.6510 


5.0798 


0.0606 


1.5382 


0.006 


6.9864 


4.7489 


0.0358 


0.9288 


0.009 


4.9630 


2.0483 


0.0778 


1.8647 



TABLE I: The parameters of the approximation Tp{k) — ap + bp tanh(dp — Cpk) for the spin system with N = 201 

depends slightly on its order in accordance with the represented formula. This conclusion is similar to that given in 
Sec lIIIAll see eq.(f20|. 

Similar to Sec lIII ATI we may introduce the cluster of MQ coherences at any time moment f as a family of such 
coherences that J^{f) > Jmin ~ 0.005. Evolution of the cluster size Nc{f) is shown in Figl5]for different N and p. 
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FIG. 5: The evolution of the cluster size Nc for iV = 201 (sohd hue), iV = 401 (dashed line), iV = 601 (dotted line) and different 
p; the widths of the lines increase with the decrease of the parameter p. The inset shows the periods of the coherence cluster 
growth and decay at small times. 

We see that the results of our simulations agree qualitatively with the experimental ones obtained in Q]. Namely, 
there is the period of the coherence cluster growth < f < 1.5 and the period of the cluster decay, r > 1.5. FiglS] 
demonstrates that the cluster size gets its maximal value at the time moment f w 1.5, which is slightly dependent on 
the both parameters N and p. This confirms our assumptions that all spins become embedded in the cluster during 
the time interval r ^ l/D, or f 1. This feature of the cluster growth in the system of equivalent spins is different 
from that in solids Fig H] demonstrates also that the maximal size of the cluster increases with the increase in N 
and slightly decreases with the increase in p. The rate of the cluster decay increases with the increase in both N and 
P- 

Comparison of Figsl3]and[5]shows that the case of the modified preparation period is more flexible in preparation 
of the coherence clusters with the desirable size because of the parameter p which does not appear in Sec lIII Al 



IV. THE CONSERVATION LAW IN THE MODEL OF THE DIPOLAR RELAXATION OF MQ NMR 

COHERENCES 

It is worth to emphasize that the appearance of MQ NMR coherences and their relaxation are determined by the 
same DDI, which is valid in the models both suggested in [ll-Q| and considered in the previous sections. This leads to 
some peculiarities of the relaxation process. We show that the sum of areas of the signals of MQ NMR coherences in 
the frequency domain is not changed in the relaxation process although their maximal amplitudes decrease. For the 
sake of simplicity, we turn to the case considered in Sec lIII Al where the decay occurs on the evolution period. Then 
the intensities of MQ NMR coherences are determined by Eq. . Performing the Fourier transform of the intensities 
Jfc(r, t) of Eq.lfTS]) over the time t of the evolution period (we suppose that Jk{T, t) = for t < and t > T, where T 
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1 '■^ 



is the duration of the evolution period) 

Jk{T,oo)^^ I ,h(T,t)e-'^'dt, (27) 
one can find that the area Afe(r) under J'k{T,uj) in the frequency domain is 

Mr)^:^ I Jk{T,u)du^ (28) 
1 1"°° 

— du: MT,t)e'''''dt = 
27r J_oo Jo 

■Mr,t)dtJ e-'^'du^-Jk{T,Q). 

Then the sum of the areas Ak for all MQ NMR coherences can be expressed as follows: 

^Afc(r) = i^jfe(r,0). (29) 

However, it is known that J/j(r, 0) = 1 [2^. Thus, ea. (p9)) means that the areas ^fc(r) arc redistributing during 
the relaxation process so that their sum is conserved. 

Similarly, replacing t with p and T with 1 in Eqs.ipT)) and (pS)) one derives the same conservation law for MQ NMR 
experiment of Sec lIII Bl 

The results of this section demonstrate some peculiarities of the used relaxation model. 



V. CONCLUSIONS 

Using the numerical methods describing the spin dynamics in large systems of equivalent spins , we study the 

decay of MQ NMR coherences in such systems. This decay is caused by the Hamiltonian Hdz appearing cither on the 
preparation or evolution period of the MQ NMR experiment. Numerical simulations are performed for the systems 
consisting of 200-600 spins. It is found that the relaxation rate of MQ NMR coherences from the highly correlated 
spin states increases with the increase in both the MQ NMR order and the number of spins. The dependence of 
the relaxation time of MQ NMR coherences on the perturbation strength p, appearing on the preparation period, is 



also investigated. We emphasize that the used model [1H3| is the first one for the experimental investigation of the 
relaxation of the correlated spin clusters of the large size. 

It is worth to note that the evolution of the intensities of MQ NMR coherences in the system of equivalent spins 
is accompanied by the reversion phenomena. Such phenomena were studied both experimentally and numerically 
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22| and the decoherence was considered as the decay of the Loschmidt echo. The reversion phenomena are not 



considered in this paper. 

All numerical simulations have been performed using the resources of the Joint Supercomputer Center (JSCC) of 
the Russian Academy of Sciences. Authors thank the anonymous referee for the valuable remarks. The work was 
supported by the Program of the Presidium of Russian Academy of Sciences No. 21 " Foundations of fundamental 
investigations of nanotechnologies and nanomaterials" . 
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